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Numerical  modeling  of  the  expansion  of  electric  thruster  plumes  provides  direct  means 
for  predicting  spacecraft  surface  contamination  and  erosion  due  to  plume  ions.  A  software 
package  named  COLISEUM  that  is  capable  of  self-consistently  modeling  plasma 
propagation  and  interactions  with  arbitrary  3-D  surfaces  is  being  developed  by  a  national 
team  of  researchers.  Despite  much  research  and  development  in  modeling  plume  expansion, 
it  is  necessary  to  continuously  validate  these  codes  using  laboratory  based  experimental  data. 
It  is  well-established  that  vacuum  chamber  facilities  affect  the  plume  of  these  devices.  Thus, 
the  models  must  not  only  describe  the  plume  expansion,  but  also  effects  of  the  vacuum 
chamber.  COLISEUM  has  been  designed  to  simulate  both  vacuum  chamber  configurations 
and  spacecraft  geometries.  This  work  presents  a  study  that  compares  results  from  a  hybrid 
particle-in-cell  model  (AQUILA)  with  Monte  Carlo  collisions  to  data  obtained  from  the 
plume  of  Busek  600 W  Hall  thruster  (BHT-HD-600).  This  data  includes  current  density,  ion 
velocities,  and  energy  distribution  data.  Also  contained  in  this  work  is  a  source  derivation 
description  from  laser  induced  fluorescence  (LIF)  data. 


I.  Introduction 

Numerical  modeling  of  the  thruster  and  surrounding  environment  provides  direct  means  for  predicting  plume 
properties  where  experimental  methods  are  limited,  such  as  predictions  of  spacecraft  interactions.  Insight  into 
the  plume  properties  and  corresponding  spacecraft  interaction  would  provide  the  community  with  a  useful  tool.  The 
Air  Force  Research  Laboratory  (AFRL)  is  leading  the  development  of  COLISEUM1,  a  3D  plasma  interaction 
framework  that  incorporates  a  plasma  expansion  tool  with  surface  interactions.  COLISEUM  has  been  designed  to 
be  usable,  flexible,  and  expandable.  It  is  capable  of  modeling  chamber  effects,  which  are  known  to  affect  the  plume 
expansion2,  as  well  as  open  boundary  conditions.  COLISEUM  has  available  any  of  four  plasma  simulation 
modules,  RAY,  PRESCRIBED_PLUME,  DRACO,  and  AQUILA.  RAY  uses  a  ray-tracing  method  to  project  a  flux 
from  a  point  surface.  PRES  CRIB  ED_PLUME  imports  and  superimposes  a  plume  distribution  onto  surfaces. 
DRACO  tracks  particles  along  a  structured  Cartesian  mesh.  AQUILA,  the  focus  of  this  study,  is  a  hybrid  particle- 
in-cell  (PIC)  model  that  tracks  particles  along  an  unstructured  tetrahedral  mesh. 
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Despite  much  research  and  development  in  modeling  plume  expansion,3’4,5  6  it  is  necessary  to  continuously 
validate  these  codes  using  laboratory  based  experimental  data.  This  paper  presents  a  study  that  compares  results 
from  an  AQUILA  simulation  to  experimental  data  from  the  BHT-HD-600  Hall  thruster.  The  numerical  study  uses  a 
simplified  Hall  thruster  to  simulate  an  electric  plume  in  a  chamber  environment.  Laser  induced  fluorescence  (LIF) 
data  is  used  to  construct  source  input  files.  Probes  in  AQUILA  collect  current  density,  ion  velocities,  and  energy 
distributions  which  are  compared  to  experimental  data  for  code  validation.  The  full  capabilities  of  AQUILA  are 
demonstrated  using  probes  to  collect  additional  information  about  the  plume. 

II.  AQUILA  code 

AQUILA,  a  hybrid  particle-in-cell  model,  has  been  developed  in  the  framework  of  COLISEUM.  AQUILA  uses 
an  unstructured  tetrahedral  mesh  to  define  surface  and  volume  geometries.  The  geometry  and  mesh  used  by 
COLISEUM  can  be  produced  using  available  commercial  modeling  and  meshing  packages.  The  mesh  can  be 
loaded  into  COLISEUM  using  any  number  of  standard  forms  including  ANSYS  and  ABAQUS.  Individual  surfaces 
are  specified  to  distinguish  surface  properties  and  define  particle/surface  interactions. 

AQUILA  contains  two  types  of  potential  solvers,  a  quasineutral  solver  and  a  non-neutral  solver.  In  this  work, 
the  quasineutral  solver  was  used.  Following  a  quasi-neutral  assumption,  the  potential  can  be  calculated  by  using  the 
inverted  Boltzmann  equation: 


kT 
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where  neo  is  a  reference  plasma  density,  (/)0  is  a  reference  electrostatic  potential,  and  Te  is  the  reference  electron 

temperature.  For  plumes  where  the  quasi-neutral  assumption  does  not  hold,  such  as  behind  a  plume  shield, 
AQUILA  contains  a  non-neutral  solver. 


Collisions  are  performed  in  AQUILA7  using  a  Direct  Simulation  Monte  Carlo  method8.  Collision  cross  sections 
for  elastic  collisions9  between  neutrals  and  ions  are  calculated  using: 
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where  vrei  is  the  relative  velocity  between  the  two  particles.  Charge  exchange  collision  cross  sections  between 
ions  and  neutrals10, 11  are  defined  as: 


crc™Xe+  =  2.2415x1  (T18  -2.766 1x10  19  log(vre/) 
= 1x10  20 (35-006  - 2-7038 iog(vre/ »2 


The  current  density  probe  in  AQUILA  tracks  ion  flux  across  a  hemisphere  at  a  user  defined  location.  The  results 
along  the  hemisphere  are  averaged  to  give  values  from  0  to  90°. 

To  decrease  computational  time,  acceleration  techniques  have  been  introduced  into  AQUILA.  One  of  the 
acceleration  techniques  utilizes  a  subcycle  routine  that  decouples  the  ion  and  neutral  movements.  Subcycling  cycles 
the  fast  particles  for  a  specified  time  step  before  moving  and  injecting  the  slow  particles  and  performing  collisions. 
Subcycling  in  AQUILA  is  shown  to  decrease  the  computational  time  for  the  simulation  to  reach  steady  state1. 
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However,  it  is  also  shown  that  too  many  subcycles  (-100)  allow  an  ion  to  leave  the  main  region  of  the  plume 
without  undergoing  any  collisions12,  an  event  that  may  not  be  entirely  physical. 

III.  Source  Definition 

Plasma  modeling  within  COLISEUM  begins  with  source  definition.  Sources  are  used  to  introduce  particles  into 
the  simulation  domain.  COLISEUM  provides  several  different  source  models,  all  specifying  a  velocity  distribution 
function  (VDF)  and  a  mass  flow  rate.  While  several  source  models  are  available,  the  test  case  used  in  this  paper 
used  the  source  FLUX_R_VZ_VR13.  As  the  name  suggests,  the  FLUX_R_VZ_VR  model  represents  the  exit  plane 
of  a  Hall  thruster  in  terms  of: 

1 .  particle  flux  versus  radial  position,  r 

2.  axial  velocity,  vz,  versus  radial  position,  r 

3.  radial  velocity,  vr,  versus  radial  position,  r 

Before  starting  the  simulation  run,  the  flux  function  is  converted  into  a  cumulative  distribution  function  for  radial 
position.  During  injection  of  particles,  the  simulation  uses  the  thruster’s  mass  flow  rate  to  determine  the  number  of 
particles  to  create  at  each  time  step.  The  previously-computed  cumulative  distribution  function  (CDF)  is  then  used  to 
place  source  particles  at  radial  distances  with  the  correct  probability.  Once  the  radial  location  of  the  injected  particle 
is  known,  the  initial  velocity  components  are  calculated  from  the  velocity  functions.  In  addition,  thermal 
components  are  added  to  these  velocities  based  on  user- specified  temperatures. 

Laser-induced  fluorescence  data  provides  a  natural  way  of  determining  the  VDF  required  by  the  simulation 
source  model.  After  processing,  the  data  taken  at  each  location  gives  a  probability  distribution  function  of  the 
velocity  component  aligned  with  the  laser  used  to  probe  the  plasma.  To  get  velocity  distribution  functions  along  a 
different  axis,  the  orientation  of  the  laser  is  changed.  For  this  source  model,  all  VDF  inputs  were  taken  from  LIF 
data14  at  an  axial  distance  of  15  mm  downstream  from  the  thruster  exit  plane.  Figure  1  shows  the  geometry  of  the 
BHT-HD-600. 


Figure  1:  BHT-HD-600  geometry  for  LIF  analysis.  The  data  for  the  LIF  analysis  is  taken  at  y  =  15  mm  which 

is  between  the  exit  plane  of  the  thruster  and  the  cathode. 

A.  Velocity  Distributions 

Figure  2  shows  a  representative  axial  VDF  for  the  BHT-HD-600  taken  at  a  radial  distance  of  24  mm  from  the 
centerline.  A  distinct  velocity  peak  is  seen  around  20,000  m/s.  To  quantify  the  mean  and  spread  of  the  peak,  the 
velocity  distribution  function  is  converted  to  a  histogram  and  a  Matlab  function  that  fits  a  specified  number  of 
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Gaussians  is  used.  Figure  3  illustrates  this  process  -  the  red  lines  show  the  fitted  Gaussians.  In  most  cases,  two 
Gaussians  were  used  to  fit  the  histogram  -  one  captures  the  property  of  the  peak  of  interest  while  the  other 
represents  background  noise  in  the  signal. 


Figure  2:  Axial  velocity  distribution  function  at  r  =  24mm  of 
BHT-HD-600. 


Figure  3:  Histogram  with  fitted  Gaussians  for  axial  velocity 
distribution  function  at  r  =  24  mm. 


A  similar  procedure  is  followed  at  each  radial  location  with  LIF  data  for  both  the  axial  and  radial  velocity 
components.  The  mean  of  the  main  Gaussian  is  used  to  represent  the  velocity  magnitude,  while  the  standard 
deviation  is  used  to  represent  the  temperature.  Since  the  source  model  assumes  azimuthally  symmetry,  only  values 
from  one  side  of  the  thruster  are  needed  -  for  this  case,  values  from  the  non-cathode  side  of  the  thruster  are  taken. 
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Upon  further  inspection  of  the  LIF  data,  it  is  observed  that  evidence  of  a  second  ion  population  is  seen  in  the 
axial  velocity  distribution  functions  found  on  the  cathode  side  of  the  thruster.  Compared  to  the  above  figures,  Figure 
4  shows  the  corresponding  axial  distribution  at  r  =  24  mm  on  the  cathode  side. 


Figure  4:  Axial  velocity  distribution  function  at  r  =  -24mm  (cathode  side)  of  BHT-HD-600. 

Two  defined  axial  velocity  peaks  are  observed,  but  when  the  radial  distributions  are  investigated,  no  difference 
in  behavior  is  seen.  It  is  conjectured  that  the  second  axial  peak  on  the  cathode  side  may  be  related  to  additional 
ionization  that  occurs  near  the  exit  plane.  Electrons  streaming  from  the  cathode  are  assumed  uniformly  azimuthally 
distributed  by  the  magnetic  field  once  they  enter  the  thruster  acceleration  channel  and  as  a  result,  ionization  is  also 
uniform.  However,  neutrals  outside  the  thruster  on  the  cathode  side  are  more  susceptible  to  ionization  by  electrons 
before  they  reach  the  channel.  Since  these  ions  are  formed  outside  of  the  main  accelerating  potential,  their  existence 
is  manifested  by  a  smaller  lower-energy  peak  in  the  axial  distributions.  As  currently  written,  the  AQUILA  source 
model  cannot  model  an  asymmetric  exit  plane  distribution.  Nevertheless,  source  model  information  including  this 
second  population  is  also  processed  in  case  one  considers  the  second  low-energy  peak  important. 

B.  Flux  Distribution 

Flux  measurements  for  the  BHT-HD-600  have  not  yet  been  completed.  An  attempt  was  made  to  extract  ion  flux 
data  from  LIF  signal  strength.  Signal  strength  in  the  linear  regime  is  proportional  to  both  ion  density  at  the  probed 
state  and  to  the  intensity  of  the  beam15.  This  approximation  holds  if  the  ion  and  electron  temperatures  and  the 
electron  density  are  relatively  constant  throughout  the  region.  The  peak  signal  strength  at  each  data  location  for  z  = 
15  mm  is  plotted  as  single  points  in  Figure  5.  Using  a  curve  fitting  algorithm,  a  Gaussian  profile  was  fit  to  this  data. 
The  flux  profile  shifts  the  center  of  the  distribution  away  from  the  centerline  of  the  channel  (at  r  =  -28  mm)  towards 
the  center  line  of  the  thruster  (r  =  24  mm).  A  flux  distribution  input  file  was  generated  using  this  distribution. 
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Figure  5:  Flux  profile  generated  from  LIF  signal  strength  data  fitted  with  a  Gaussian  distribution.  The 
Gaussian  fit  to  the  LIF  data  is  centered  at  r  =  -24  mm  rather  than  over  the  centerline  of  the  thruster  at  r  =  -28 

mm. 


IV.  Problem  Description 

The  chamber  and  thruster  geometry  for  the  simulations  is  based  on  the  BHT-HD-600  Hall  thruster  live  tests 
inside  Chamber  6  at  the  Air  Force  Research  Laboratory.  During  the  thruster  testing,  graphite  panels  were  added  to 
the  chamber  to  lower  sputter  rates  of  the  chamber  and  reduce  re-deposition  of  the  chamber  materials  back  onto  the 
thruster.  The  simulation  preserved  the  placement  and  orientation  of  the  graphite  panels  while  simplifying  the 
geometry  of  the  panels  and  other  components  to  increase  mesh  quality.  Figure  6  shows  the  surface  mesh  of  the 
chamber,  graphite  panels,  and  thruster  orientation.  The  thruster  fires  in  the  negative  z  direction. 

The  chamber  modeled  is  a  1.8  m  diameter,  2.9  m  length  stainless  steel  chamber.  The  graphite  panels  are 
modeled  as  semi-circular  shapes  concentric  with  the  chamber,  60  cm  in  width  and  lengths  from  79  cm  to  1.22  m. 
Cryogenic  pumps  are  also  included  in  the  geometry  at  the  rear  of  the  chamber,  87  cm  wide  and  83  cm  apart.  The 
chamber  background  pressure  was  set  to  7  x  10'4  Pa,  corrected  for  Xe. 
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z  (m) 


Figure  6:  Chamber  6  geometry  and  mesh.  The  green  elements  are  the  graphite  panels.  The  blue  elements 
are  the  cryogenic  pumps.  The  thruster  is  the  orange  box  in  the  center  of  the  chamber. 


The  simulated  Hall  thruster  geometry,  shown  in  Figure  7,  is  the  plasma  emitting  source  in  the  simulation. 
Particles  are  emitted  from  the  annular  region  of  the  thruster  face.  Based  on  the  geometry  of  the  BHT-HD-600,  the 
annulus  has  an  inner  diameter  of  24  mm  and  an  outer  diameter  of  32  mm. 


Figure  7:  BHT-HD-600  simplified  mesh.  The  yellow  elements  are  defined  as  the  particle  emitting  source. 

Notice  the  change  in  element  resolution  from  the  front  face  to  the  back. 

The  volume  mesh  used  in  this  simulation  is  shown  in  Figure  8.  The  mesh  is  a  structured  mesh  in  the  core  region 
of  the  plume,  produced  to  increase/control  the  resolution  of  elements  in  the  immediate  region  of  the  plume.  This 
region  is  defined  by  the  annulus  extending  outward  30  cm  from  the  thruster  face,  with  the  end  elements  triple  the 
size  of  the  original  elements.  An  unstructured  mesh  fills  in  the  remainder  of  the  volume. 
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Figure  8:  The  simulated  volume  mesh  with  a  structured  mesh  in  the  region  of  the  plume  and  an 
unstructured  mesh  in  the  remaining  volume.  The  structured  mesh  is  used  to  control  the  resolution  of  the 

mesh  in  the  immediate  region  of  the  plume. 


The  annular  region,  shown  in  Figure  7,  is  the  particle  emitting  source  of  the  simulation.  Based  on  experimental 
findings16,  the  particles  emitted  include  not  only  Xe  neutrals  and  singly  charged  Xe  ions,  but  multiply  charged  Xe 
ions  as  well.  The  distribution  of  neutrals  is  defined  by  a  drifting  Maxwellian  source  model  with  a  thermal  drift 
velocity  of  297  m/s  and  a  temperature  of  700  K.  The  mass  flow  rate  of  neutral  particles  into  the  simulation  reflects 
10%  of  the  unionized  propellant.  Source  files  for  doubly  charged  ions  were  generated  using  the  singly  charged  ion 
source  files,  increasing  the  magnitude  of  velocity  by  the  V 2  to  represent  an  increase  in  speed  related  to  their  charge 
according  to  the  energy  relation: 


-mv2  =  Zqkcj) 


where  Z  is  the  ion  charge,  q  is  the  elementary  charge,  and  (/)  is  the  potential. 


Although  evidence  of  a  low  velocity  ion  population  was  seen  in  the  axial  LIF  data,  this  population  was  only 
evidenced  in  the  cathode  region  of  the  plume.  Since  the  majority  of  the  plume  lacks  this  second  population,  the 
simulation  was  completed  using  only  the  main  population  of  ions,  and  the  flux  profile  shown  in  Figure  5.  Ion 
temperatures  of  0.81  eV  in  the  axial  and  azimuthal  directions  and  0.81  eV  in  the  radial  direction  were  used. 


The  electrostatic  potential  is  computed  using  the  quasi-neutral  approach  in  AQUILA.  The  reference 
potential^ is  set  to  40  eV,  and  the  reference  density  neo  is  sampled  from  the  domain  at  a  point  centered  on  the 
annular  channel.  The  electron  temperature  T  is  specified  using  a  polytropic  model: 


T  =  T 


\UeoJ 


r- 1 


where  the  reference  electron  temperature  Te  is  set  to  10.0  eV  and  y  is  1.3. 


The  timing  scheme  in  AQUILA  is  used  to  determine  when  the  simulation  reaches  a  steady  state.  Subcycling  was 
implemented  in  this  simulation  to  decrease  computation  time.  Ions  were  subcycled  10  times  (with  an  ion  time  step 
of  2  x  10-7  s)  every  cycle.  The  number  of  cycles  completed  was  2500.  The  region  of  the  domain  specified  as  the 
cryogenic  pumps  takes  incoming  particles  out  of  the  domain  at  a  user  specified  rate.  In  this  model,  a  cryogenic 
sticking  coefficient  of  30%  was  chosen,  meaning  that  30%  of  particles  hitting  the  pumps  will  be  removed  from  the 
simulation. 
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V.  Comparison  against  Experimental  Data 

The  plume  studies  of  the  BHT-HD-600  conducted  by  Ekholm,  et  al16,  provided  measurements  of  the  ion  current 
density  profile,  ion  energy  distributions,  and  ion  species  fraction  distributions  using  a  nude  Faraday  probe,  retarding 
potential  analyzer  (RPA),  and  ExB  probe.  LIF  measurements  were  taken  from  Charles  and  Hargus14.  This  suite  of 
data  serves  as  a  comparison  to  the  test  cases  completed  using  AQUILA  and  the  input  parameters  defined  above. 

A.  Current  Density 

The  current  density  probe  in  AQUILA  was  set  up  to  sample  particles  along  a  hemisphere  60  cm  in  front  of  the 
thruster,  allocating  200  bins  for  sampling.  Figure  9  shows  the  current  density  from  the  test  case  compared  to 
experimental  results  from  Faraday  probe  measurements  at  60  cm.  As  shown  in  Figure  9,  the  trend  in  the  current 
density  profile  agrees  with  experimental  data  across  the  region,  although  the  magnitudes  only  agree  from  -40°  to  40°. 
In  the  wing  region  of  the  plume,  the  data  for  the  test  case  is  lower  than  what  is  seen  experimentally.  These  results 
indicate  a  narrow  beam  divergence.  Experiments  are  known  to  have  difficulties  in  this  region  due  to  charge 
exchange  collisions  and  secondary  electron  emission  from  ion  impact  with  the  probe.  The  difference  could  also 
indicate  greater  collision  rates  in  the  plume  than  what  is  modeled.  Increasing  collision  rates  would  lead  to  greater 
plume  divergence  and  spread  more  ions  out  towards  the  wings.  Despite  the  differences,  the  similar  trend  between 
the  two  sets  of  data  suggests  that  only  minor  refinements  are  needed  in  the  model  to  bring  the  both  sets  of  data  into 
agreement. 


Figure  9:  Current  density  at  60  cm  from  Ekholm,  et  al  and  from  AQUILA  test  case.  The  trends  in  both 
sets  of  data  are  similar,  with  the  greatest  differences  occurring  in  the  wings  of  the  plume  from  +40°  to  +90°. 
Greater  collision  rates  in  the  model  would  broaden  the  spread  of  the  test  case  profile. 


B.  Near  Field  Ion  Velocities 

Figure  10  shows  a  comparison  of  axial  ion  velocities  at  x  =  0  mm,  z  =  60  mm  from  the  thruster  exit  with  LIF 
peak  data  values  taken  by  Charles  and  Hargus  and  an  uncertainty  of  500  m/s.  Figure  11  is  the  radial  velocity 
comparison  in  this  same  region.  Both  figures  depict  actual  particle  velocities  for  the  model  data,  similar  to  a 
snapshot  taken  of  the  flow.  AQUILA  tracks  particles  based  on  the  collisions  they  undergo  and  attaches  a  prefix  to 
the  species  name.  An  SRC  delimiter  signifies  a  beam  ion  that  has  not  experienced  any  collisions.  EL_COL  and 
CEX_COL  label  ions  that  have  undergone  elastic  and  charge  exchange  collisions,  respectively.  It  can  be  noted  that 
the  trend  of  the  simulated  source  ions  follows  the  trend  of  the  LIF  data,  with  the  axial  velocities  around  20  km/s  and 
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a  3  km/s  deep  depression  around  y  =  0  mm.  The  doubly  charged  source  ions  (indicated  as  SRC_XE++)  have 
velocities  around  28,000  km/s,  as  they  were  emitted  from  the  source  at  a  greater  velocity  than  their  singly  charged 
counterparts.  There  is  a  large  concentration  of  source  ions  at  y  =  0  mm,  indicating  a  significant  number  of  ions  from 
opposite  sides  of  the  channel  crossing  over  the  centerline.  This  same  trend  is  also  seen  in  the  radial  data.  The  radial 
velocities  follow  a  linear  pattern  across  both  sections  of  the  channel,  changing  from  a  divergent  to  convergent 
profile  at  y  =  +-30  mm,  noted  by  a  shift  in  velocities  from  positive  to  negative.  As  is  expected,  the  ions  undergoing 
collisions  are  shown  to  have  axial  velocities  around  1000  m/s  and  radial  velocities  near  0  m/s.  There  is  a  significant 
amount  of  spread  in  the  velocity  distribution  which  is  also  seen  in  data  taken  by  Charles  and  Hargus.  While  the 
greatest  population  of  ions  occurs  at  the  peak,  higher  velocity  and  lower  velocity  populations  exist.  Figures  10  and 
1 1  suggest  that  the  model  portrays  consistent  results  in  the  near-field  as  compared  to  experimental  data. 
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Figure  10:  Axial  LIF  comparison  at  60  mm  in  front  of  thruster.  The  trends  from  both  the  experimental 
and  test  case  data  agree,  with  a  significant  amount  of  distribution  in  the  test  case  data  also  seen  in  the  data  by 

Charles  and  Hargus. 
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at  x=  0  mm,  z  =  60  mm 


Figure  11:  Radial  LIF  comparison  at  60  mm  in  front  of  the  thruster.  The  trends  from  both  the 
experimental  and  test  case  data  agree.  Note  the  change  of  ion  velocities  from  convergent  to  divergent  flow  at 
y  =-+30  mm,  when  they  cross  from  negative  to  positive  velocities. 

C.  Far  Field  Energy  Distribution 

The  energy  distribution  of  ions  in  the  simulation  can  be  determined  from  the  ion  velocities.  Using  the  equation 
for  kinetic  energy:  KE  =  Vi  mv2,  with  the  mass  of  a  Xe  ion  as  2.182  xlO'25  kg,  the  kinetic  energy  per  charge  is 
computed  along  a  hemisphere  60  cm  in  front  of  the  thruster,  at  specific  angles  (0°,  ±30°,  ±60°,  ±90°)  from  the 
thruster  centerline  (0°).  The  resulting  plots,  normalized  by  area,  are  shown  in  Figure  12-15  compared  against  RPA 
experimental  data.  The  RPA  probe  measures  an  energy  distribution  per  charge;  the  peak  values  shown  in  Figures 
12-15  are  the  most  probable  ion  energy  levels.  Along  the  centerline  of  the  thruster,  the  most  probable  energy  level 
the  model  predicts  is  277  eV/q,  compared  to  the  slightly  lower  experimental  measurement  of  262  eV/q.  Figure  12 
also  shows  that  the  model  does  not  display  a  significant  number  of  low  energy  ions  around  the  centerline,  as  can  be 
noticed  in  the  data  by  Ekholm,  et  al.  Figure  13  displays  the  energy  distribution  at  ±30°.  The  AQUIFA  test  case 
shows  a  slightly  lower  high  energy  peak,  at  248  eV/q  compared  to  the  experimental  data,  which  is  278  eV/q.  The 
presence  of  low  energy  ions  (around  34  eV/q)  begins  to  become  evident  at  ±30°,  although  not  as  prominent  as  the 
experimental  RPA  data.  At  ±60°,  the  low  energy  ions  dominate  the  AQUIFA  test  case  energy  distribution,  as  seen 
in  Figure  14;  no  evidence  of  any  high  energy  ions  is  seen.  Although  a  slight  peak  is  seen  around  160  eV/q  for  the 
model  data  at  +60°,  the  great  majority  of  the  ion  energy  remains  around  30  eV/q.  This  trend  differs  significantly 
from  the  experimental  data,  where  the  energy  distribution  is  divided  between  high  energy  (282  eV/q)  and  low 
energy  (51  eV/q)  particles.  Figure  15  shows  the  ion  energy  distribution  in  the  wings  of  the  plume  at  ±90°.  The 
AQUIFA  test  case  data  declines  earlier  than  the  experimental  data,  but  both  sets  agree  that  the  majority  of  the  ions 
have  an  energy  per  charge  of  36  eV/q.  The  graphs  in  Figures  13-15  are  plotted  to  show  symmetry  across  the 
centerline  of  the  thruster.  While  minor  discrepancies  exist,  the  data  shows  good  agreement  between  opposite  sides 
of  the  thruster. 

In  a  Hall  thruster  plume,  the  high  energy  ions  are  typically  beam  ions  while  the  low  energy  ions  are  formed  from 
charge  exchange  collisions.  The  ions  residing  in  the  middle  region  are  usually  due  to  elastic  scatter  collisions.  The 
experimental  data  shows  evidence  of  beam  ions  at  ±60°.  This  is  not  seen  in  the  model,  suggesting  that  the  model 
depicts  a  narrow  beam  divergence.  The  current  density  plot  (shown  in  Figure  9)  also  suggests  this.  Beam 
divergence  in  the  model  is  controlled  by  four  parameters:  a)  the  radial  velocity  source  files  b)  charge  exchange 
collisions  and  to  a  lesser  extent  elastic  collisions  c)  source  ion  temperature  and  iv)  the  self-consistent  E-field.  From 
the  FIF  comparison  at  60  mm,  in  Figures  10  and  11,  it  can  be  inferred  that  the  source  model  is  satisfactory.  This 
indicates  that  the  two  remaining  parameters,  the  collision  model  and  the  E-field,  may  be  the  principle  contributors  to 
the  narrow  beam  divergence.  The  absence  of  low  energy  ions  (created  by  collisions)  in  much  of  the  model  data 
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supports  this  theory.  Additional  investigation  into  the  AQUILA  collision  model  and  the  potential  solver  may  bring 
both  sets  of  data  into  agreement. 


Energy  per  Charge  (eV/q) 

Figure  12:  Ion  energy  distribution  at  r  =  60  cm  from  Ekholm,  et  al  and  from  the  AQUILA  test  case,  along 
the  centerline  of  the  thruster.  The  peak  energy  value  for  the  model  (shown  in  blue)  is  272  eV/q,  while  the 

RPA  data  has  a  peak  value  of  262  eV/q. 
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Figure  13:  Ion  energy  distribution  at  r  =  60  cm  from  Ekholm,  et  al  and  from  the  AQUILA  test  case,  at 
+30°.  The  peak  energy  values  for  the  model  (248  eV/q)  are  slightly  lower  than  the  RPA  data  (278  eV/q).  Also, 
the  absence  of  low  energy  ions  can  be  noted  between  0  eV/q  and  200  eV/q. 


Figure  14:  Ion  energy  distribution  at  r  =  60  cm  from  Ekholm,  et  al  and  from  the  AQUILA  test  case,  at 
+60°.  The  model  shows  no  middle  to  high  energy  level  ions,  unlike  the  RPA  data.  The  peak  energy  value  in 
the  model  is  34  eV/q,  where  the  RPA  data  shows  peak  distributions  at  282  eV/q  and  50  eV/q. 
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Energy  per  Charge  (eV/q) 


Figure  15:  Ion  energy  distribution  at  r  =  60  cm  from  Ekholm,  et  al  and  from  the  AQUILA  test  case  at 
±90°.  The  model  and  RPA  data  agree  with  a  peak  energy  of  36  eV/q,  although  the  model  data  declines  sooner 

than  the  experimental  data. 


VI.  Results 

Computer  models  provide  an  insight  into  the  less  tangible  characteristics  of  the  plume.  Using  the  test  case 
above,  plots  were  constructed  for  several  properties  including  plasma  potential,  electron  temperature,  plasma 
density,  and  Debye  length.  These  plots  are  shown  in  Figures  16-23  in  both  a  full  chamber  perspective  and  a  close  up 
view  of  the  thruster.  Potential  plots  are  shown  in  Figures  16  and  17.  As  is  predicted,  the  greatest  potential  (48  V)  is 
seen  directly  in  front  of  the  thruster,  expanding  radially  outwards  and  decreasing  with  increasing  distance  from  the 
thruster.  In  the  detail  view,  the  potential  contours  converge  along  the  thruster  centerline  a  distance  out  from  the 
thruster  face.  Electron  temperature,  shown  in  Figures  18  and  19  follow  the  same  trend,  with  a  maximum 
temperature  of  10  eV  directly  in  front  of  the  thruster.  Plasma  density,  illustrated  in  Figures  20  and  21  show  the 
highest  concentration  of  ions  (3.6  x  1017/m3)  directly  in  front  of  the  thruster,  dropping  off  rapidly  outside  the 
immediate  region  of  the  plume.  Debye  length  is  shown  in  Figures  22  and  23.  The  Debye  length  inside  the  chamber 
rises  from  6.5  x  10'5m  in  front  of  the  thruster  to  1.2  x  10'3  m  at  the  back  of  the  chamber. 
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Figure  16:  Chamber  view  of  electrostatic  potential.  The  potential  profile  diverges  radially  outward  in  the  main 
region  of  the  beam,  decreasing  in  magnitude  from  48  V  to  0  V  at  the  back  of  the  chamber. 


Figure  17:  Thruster  view  of  electrostatic  potential.  The  potential  quickly  drops  off  from  48  V  in  front  of  the 
thruster  and  converges  along  the  centerline  of  the  thruster  downstream. 
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Figure  18:  Chamber  view  of  electron  temperature.  The  temperature  ranges  from  10  eV  in  front  of  the 
thruster  to  0.5  eV  in  the  back  of  the  chamber.  The  profile  radiates  outward  following  the  same  trend  seen  in 

the  potential. 


Figure  19:  Thruster  view  of  the  electron  temperature.  The  electron  temperature  matches  the  potential  profile 
by  converging  along  the  centerline  downstream  of  the  thruster. 
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Figure  20:  Chamber  view  of  plasma  density.  The  greatest  concentration  of  plasma  is  directly  in  front  of  the 
thruster.  The  plasma  density  drops  outside  the  main  region  of  the  plume  to  a  background  density  of  2  x  1016. 


Figure  21:  Thruster  view  of  plasma  density.  The  plasma  density  also  converges  on  the  centerline  of  the 
thruster,  similar  to  the  potential  and  electron  temperature  profiles  as  seen  in  Figures  17  and  19. 
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Figure  22:  Chamber  view  of  Debye  length.  The  Debye  length  is  8.0  x  10'4  in  the  wings  of  the  thruster  and 

rises  to  1.2  x  10'3  in  the  far  reaches  of  the  chamber. 


Figure  23:  Thruster  view  of  Debye  length.  The  Debye  length  in  the  front  of  the  thruster  is  6.5  x  10'5  and  1.1  x 

10'4  downstream,  along  the  centerline  of  the  thruster. 
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VII.  Conclusion 

This  particular  test  case  has  shown  that  the  results  from  AQUILA  agree  with  experimental  data  in  magnitude  and 
trend  while  not  always  agreeing  in  distributions.  Results  from  near-field  velocity  profiles  indicate  that  the  source 
model  velocity  distributions  are  consistent  with  experimental  data.  The  comparisons  against  current  density  and 
energy  distribution  data  indicate  a  narrow  beam  divergence,  possibly  caused  by  low  collision  rates  and  an 
inconsistent  E-field.  Future  work  needs  to  be  completed  in  AQUILA  to  determine  which  parameters  effectively 
decide  beam  divergence,  starting  with  the  collision  model  and  the  potential  solver. 

VIII.  F  utur  e  W  ork 

This  work  is  presented  using  only  one,  optimal  test  case.  The  results  indicate  that  adjustments  to  the  collision 
model  and  potential  solver  may  yield  more  accurate  results.  Also,  by  changing  other  parameters  of  the  simulation 
such  as  background  pressure,  mesh  quality,  and  source  file  inputs,  the  results  from  the  simulation  could  change.  A 
sensitivity  analysis  of  AQUILA  to  input  parameters  is  of  future  interest.  In  addition,  more  data  from  the  BHT-HD- 
600  may  become  available  in  the  future.  At  that  time,  it  would  be  pertinent  to  compare  the  data  to  results  from  this 
model. 
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